/**
 * 
 */
package cn.cityhouse.avm.utils;

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics;

import cn.cityhouse.avm.stat.SummaryStats;

/**
 * This class implements some useful <b>GIS util functions</b>
 * 
 * @author <a href="mailto:wangcl@lreis.ac.cn">ChenLiang Wang</a>
 *
 */
public class GISUtils {

	/**
	 * @deprecated function computing great circle distance for coordinates
	 *             (WGS84)
	 * @param p1
	 *            double array of the first coordinates
	 * @param p2
	 *            double array of the second coordinates
	 * @return great circle distance
	 */
	@Deprecated
	public static double distance(double[] p1, double[] p2) {
		double sum = 0;
		for (int i = 0; i < p1.length; i++) {
			final double dp = p1[i] - p2[i];
			sum += dp * dp;
		}
		return Math.sqrt(sum);
	}

	/**
	 * @deprecated Calculate distance between two points in latitude and
	 *             longitude taking into account height difference. If you are
	 *             not interested in height difference pass 0.0. Uses Haversine
	 *             method as its base.
	 * @param coord_pt1
	 *            Start point
	 * @param coord_pt2
	 *            End point
	 * @param el1
	 *            Start altitude in meters
	 * @param el2
	 *            End altitude in meters
	 * @return Distance in Meters
	 */
	@Deprecated
	public static double distance(double[] coord_pt1, double[] coord_pt2, double el1, double el2) {

		final int R = 6371; // Radius of the earth
		double lon1 = coord_pt1[0];
		double lat1 = coord_pt1[1];
		double lon2 = coord_pt2[0];
		double lat2 = coord_pt2[1];
		Double latDistance = Math.toRadians(lat2 - lat1);
		Double lonDistance = Math.toRadians(lon2 - lon1);
		Double a = Math.sin(latDistance / 2) * Math.sin(latDistance / 2) + Math.cos(Math.toRadians(lat1))
				* Math.cos(Math.toRadians(lat2)) * Math.sin(lonDistance / 2) * Math.sin(lonDistance / 2);
		Double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
		double distance = R * c * 1000; // convert to meters

		double height = el1 - el2;

		distance = Math.pow(distance, 2) + Math.pow(height, 2);

		return Math.sqrt(distance);
	}

	/**
	 * function computing great circle distance for coordinates (WGS84)
	 * 
	 * @param coord_pt1
	 *            double array of the first coordinates
	 * @param coord_pt2
	 *            double array of the second coordinates
	 * @return Great Circle Distance
	 */
	public static double distance_84(double[] coord_pt1, double[] coord_pt2) {
		double F, G, L, sinG2, cosG2, sinF2, cosF2, sinL2, cosL2, S, C;
		double w, R, a, f, D, H1, H2;
		double lat1R, lat2R, lon1R, lon2R, DE2RA;
		DE2RA = Math.PI / 180;
		a = 6378.137; /* WGS-84 equatorial radius in km */
		f = 1.0 / 298.257223563; /* WGS-84 ellipsoid flattening factor */
		if (Double.valueOf(coord_pt1[1]).equals(coord_pt2[1])) {
			if (Double.valueOf(coord_pt1[0]).equals(coord_pt2[0])) {
				return 0d;
			} else if (Double.valueOf(Math.abs(Math.abs(coord_pt1[0]) + Math.abs(coord_pt2[0]))).equals(360)) {
				return 0d;
			}
		}
		lat1R = coord_pt1[1] * DE2RA;
		lat2R = coord_pt2[1] * DE2RA;
		lon1R = coord_pt1[0] * DE2RA;
		lon2R = coord_pt2[0] * DE2RA;
		F = (lat1R + lat2R) / 2.0;
		G = (lat1R - lat2R) / 2.0;
		L = (lon1R - lon2R) / 2.0;
		sinG2 = Math.pow(Math.sin(G), 2);
		cosG2 = Math.pow(Math.cos(G), 2);
		sinF2 = Math.pow(Math.sin(F), 2);
		cosF2 = Math.pow(Math.cos(F), 2);
		sinL2 = Math.pow(Math.sin(L), 2);
		cosL2 = Math.pow(Math.cos(L), 2);

		S = sinG2 * cosL2 + cosF2 * sinL2;
		C = cosG2 * cosL2 + sinF2 * sinL2;

		w = Math.atan(Math.sqrt(S / C));
		R = Math.sqrt(S * C) / w;

		D = 2 * w * a;
		H1 = (3 * R - 1) / (2 * C);
		H2 = (3 * R + 2) / (2 * S);

		return D * (1 + f * H1 * sinF2 * cosG2 - f * H2 * cosF2 * sinG2);
	}

	/**
	 * @param pts
	 *            the coordinates in the entire study area
	 * @return the rectangle coordinates of the study area
	 */
	public static double[][] getBbox(double[][] pts) {
		RealMatrix pts_coord = new Array2DRowRealMatrix(pts);
		DescriptiveStatistics descx = SummaryStats.summaryStats(pts_coord.getColumn(0), false);
		DescriptiveStatistics descy = SummaryStats.summaryStats(pts_coord.getColumn(1), false);
		double[][] d_bbox = new double[][] { { descx.getMin(), descx.getMax() }, { descy.getMin(), descy.getMax() } };
		return (d_bbox);
	}

	/**
	 * @param coords
	 *            the input coordinates
	 * @return boolean indicating whether the input coordinates couble be
	 *         geographic
	 */
	public static boolean isLonglat(RealMatrix coords) {
		if (coords == null || coords.getColumnDimension() <= 0 || coords.getRowDimension() <= 0) {
			throw new IllegalArgumentException();
		}
		DescriptiveStatistics ux = SummaryStats.summaryStats(coords.getColumn(0), false);
		DescriptiveStatistics uy = SummaryStats.summaryStats(coords.getColumn(1), false);
		return (ux.getMin() >= -180 && ux.getMax() <= 180 && uy.getMin() >= -90 && uy.getMax() <= 90);
	}

	private GISUtils() {
		// TODO Auto-generated constructor stub
		throw new Error("Do not instantiating this Class");
	}

}
